R.version
## _
## platform x86_64-apple-darwin20
## arch x86_64
## os darwin20
## system x86_64, darwin20
## status
## major 4
## minor 4.1
## year 2024
## month 06
## day 14
## svn rev 86737
## language R
## version.string R version 4.4.1 (2024-06-14)
## nickname Race for Your Life
library(knitr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(here)
## here() starts at /Users/yxy/Dropbox/Collaborations/3006_iterative
library(broom)
library(broom.mixed)
library(writexl)
library(ggpubr)
library(gghighlight)
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(sjPlot)
library(performance)
library(RColorBrewer)
color_list <- c(
"#FF9646", # orange for VH
"#5FB991" # green for VD
)
exp1_data_original <- read_csv(here("new_analysis_202406", "data_for_new_analysis", "cleaned_data", "exp1_new_cleaned_data_20240607.csv"))
## Rows: 9108 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (18): condition, participant, phase, training_item, training_stem, train...
## dbl (4): training_loop, training_trial_num, training_rt, testing_trial_num
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
exp1_data_original
exp1_data_good <- exp1_data_original %>%
mutate(
# add a column for each of the three choices: yes ~ 1, no ~ 0
harmonic_index = case_when(
response_type == "harmonic" ~ 1,
TRUE ~ 0
),
disharmonic_index = case_when(
response_type == "disharmonic" ~ 1,
TRUE ~ 0
),
unseen_index = case_when(
response_type == "unseen" ~ 1,
TRUE ~ 0
)
) %>%
group_by(participant) %>%
mutate(unseen_rate = sum(unseen_index) / 42) %>%
ungroup() %>%
filter(unseen_rate < 1/3)
exp1_data_good
exp1_data_original %>% select(condition, participant) %>%
distinct() %>%
group_by(condition) %>%
count()
exp1_data_good %>% select(condition, participant) %>%
distinct() %>%
group_by(condition) %>%
count()
exp1_unseen_only <- exp1_data_good %>%
filter(response_type == "unseen") %>%
mutate(
unseen_dominant_index = case_when(
condition == "vh" & unseen_type == "harmonic" ~ 1,
condition == "vd" & unseen_type == "disharmonic" ~ 1,
TRUE ~ 0
))
exp1_unseen_only
No difference between two conditions; the unseen choice trials can be
excluded
exp1_unseen_only %>%
group_by(condition) %>%
summarize(total = length(unseen_type),
n_dominant_unseen = sum(unseen_dominant_index),
dominant_unseen_rate = n_dominant_unseen / total)
t.test(unseen_dominant_index ~ condition, data = exp1_unseen_only) %>%
tidy()
Excluding the unseen choice trial;
Marking the dominant choices.
exp1_data <- exp1_data_good %>%
filter(response_type != "unseen") %>%
mutate(
dominant_index = case_when(
condition == "vh" & response_type == "harmonic" ~ 1,
condition == "vd" & response_type == "disharmonic" ~ 1,
TRUE ~ 0
)
)
exp1_data
Setting the reference level: VH
exp1_data$condition <- factor(exp1_data$condition, levels = c("vd", "vh"))
Setting up sum coding scheme
contrasts(exp1_data$condition) <- contr.sum(2)
contrasts(exp1_data$condition)
## [,1]
## vd 1
## vh -1
exp1_model <- glmer(data = exp1_data,
dominant_index ~ condition + (1 | participant) + (1 | testing_stem),
family = "binomial")
summary(exp1_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: dominant_index ~ condition + (1 | participant) + (1 | testing_stem)
## Data: exp1_data
##
## AIC BIC logLik deviance df.resid
## 2840.4 2863.0 -1416.2 2832.4 2123
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9012 -1.0366 0.6430 0.8287 1.3620
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant (Intercept) 0.1200 0.3464
## testing_stem (Intercept) 0.0524 0.2289
## Number of obs: 2127, groups: participant, 56; testing_stem, 42
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.37372 0.07426 5.033 4.84e-07 ***
## condition1 -0.28269 0.06529 -4.330 1.49e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## condition1 -0.096
The VH dominant condition
exp1_data_vh <- exp1_data %>% filter(condition == "vh")
t.test(exp1_data_vh$dominant_index, mu = 0.72) %>% tidy()
The VD dominant condition
exp1_data_vd <- exp1_data %>% filter(condition == "vd")
t.test(exp1_data_vd$dominant_index, mu = 0.72) %>% tidy()
exp1_summary <- emmeans(exp1_model, pairwise ~ condition, type = "response")$emmeans %>%
tidy(conf.int = TRUE)
exp1_summary
Setting the level order
exp1_summary$condition <- factor(exp1_summary$condition, levels = c("vh", "vd"))
exp1_plot <- ggplot(
data = exp1_summary,
aes(
x = condition,
y = prob
)
) +
geom_bar(
stat = "identity",
aes(
fill = condition
),
color = "black",
width = 0.8,
alpha = 0.85
) +
geom_errorbar(
aes(
ymin = conf.low,
ymax = conf.high
),
width = 0.2
) +
geom_hline(yintercept = 0.5,
linetype = "dashed",
size = 0.75,
color = "grey70") +
geom_hline(yintercept = 0.72,
linetype = "dashed",
size = 0.75,
color = "grey70") +
geom_signif(
comparisons = list(c("vh", "vd")),
y_position = 0.8,
tip_length = 0.2,
annotations = "***",
textsize = 12
) +
scale_fill_manual(values = c(color_list[1], color_list[2])) +
scale_x_discrete(limits = c("vh", "vd"),
labels = c("Harmony-\ndominant", "Disharmony-\ndominant")) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.25),
labels = scales::percent
) +
labs(x = "Condition",
y = "% Dominant pattern") +
theme_bw(base_size = 25) +
theme(
legend.position = "none",
axis.text.x = element_text(colour = "black"),
axis.text.y = element_text(colour = "black")
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
exp1_plot
ggsave(plot = exp1_plot,
filename = here("new_analysis_202406", "figures", "exp1_plot.pdf"),
height = 6,
width = 6)
ggsave(plot = exp1_plot,
filename = here("new_analysis_202406", "figures", "exp1_plot.png"),
dpi = 300,
height = 6,
width = 6)
exp2_data_original <- read_csv(here("new_analysis_202406", "data_for_new_analysis", "cleaned_data", "exp2_new_cleaned_data_20240610.csv"))
## Rows: 7866 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (12): condition, participant, phase, pre_training_item, testing_stem, pa...
## dbl (4): order, dominant_index, unseen_index, unseen_dominant_index
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
exp2_data_original
The rates in both the pre-contact testing and post-contact testing
phases should be lower than 1/3
exp2_data_good <- exp2_data_original %>%
mutate(
pre_testing_unseen_index = case_when(
phase == "pre_testing" & unseen_index == 1 ~ 1,
TRUE ~ 0
),
post_testing_unseen_index = case_when(
phase == "post_testing" & unseen_index == 1 ~ 1,
TRUE ~ 0
)
) %>%
group_by(participant) %>%
mutate(
pre_testing_unseen_rate = sum(pre_testing_unseen_index) / 8,
post_testing_unseen_rate = sum(post_testing_unseen_index) / 34
) %>%
ungroup() %>%
filter(
pre_testing_unseen_rate < 1/3 & post_testing_unseen_rate < 1/3
)
exp2_data_good
exp2_data_original %>% select(condition, participant) %>%
distinct() %>%
group_by(condition) %>%
count()
exp2_data_good %>% select(condition, participant) %>%
distinct() %>%
group_by(condition) %>%
count()
exp2_unseen_only <- exp2_data_good %>%
filter(unseen_index == 1)
exp2_unseen_only
No difference between two conditions; the unseen choice trials can be
excluded
exp2_unseen_only %>%
filter(phase == "pre_testing") %>%
group_by(condition) %>%
summarize(total = length(unseen_harmony),
n_dominant_unseen = sum(unseen_dominant_index),
dominant_unseen_rate = n_dominant_unseen / total)
t.test(unseen_dominant_index ~ condition, data = exp2_unseen_only %>% filter(phase == "pre_testing")) %>%
tidy()
No difference between two conditions; the unseen choice trials can be
excluded
exp2_unseen_only %>%
filter(phase == "post_testing") %>%
group_by(condition) %>%
summarize(total = length(unseen_harmony),
n_dominant_unseen = sum(unseen_dominant_index),
dominant_unseen_rate = n_dominant_unseen / total)
t.test(unseen_dominant_index ~ condition, data = exp2_unseen_only %>% filter(phase == "post_testing")) %>%
tidy()
exp2_data_contact <- exp2_data_good %>%
filter(phase == "contact")
exp2_data_contact
Excluding the unseen choice trial
exp2_data_testing <- exp2_data_good %>%
filter(phase != "pre_training") %>%
filter(phase != "contact") %>%
filter(unseen_index == 0)
exp2_data_testing
Setting the reference level: VH
exp2_data_contact$condition <- factor(exp2_data_contact$condition, levels = c("vd", "vh"))
Setting up sum coding scheme
contrasts(exp2_data_contact$condition) <- contr.sum(2)
contrasts(exp2_data_contact$condition)
## [,1]
## vd 1
## vh -1
exp2_model_contact <- glmer(data = exp2_data_contact,
dominant_index ~ order * condition
+ (1 | participant) + (0 + order | participant) + (1 | testing_stem),
family = "binomial")
summary(exp2_model_contact)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: dominant_index ~ order * condition + (1 | participant) + (0 +
## order | participant) + (1 | testing_stem)
## Data: exp2_data_contact
##
## AIC BIC logLik deviance df.resid
## 2920.4 2960.3 -1453.2 2906.4 2201
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6341 -1.0796 0.6014 0.8551 1.1590
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant (Intercept) 0.0053188 0.07293
## participant.1 order 0.0001988 0.01410
## testing_stem (Intercept) 0.0384349 0.19605
## Number of obs: 2208, groups: participant, 46; testing_stem, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.066390 0.102054 0.651 0.515
## order 0.015239 0.003893 3.915 9.04e-05 ***
## condition1 -0.012193 0.089312 -0.137 0.891
## order:condition1 -0.006116 0.003869 -1.581 0.114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) order cndtn1
## order -0.634
## condition1 -0.017 0.025
## ordr:cndtn1 0.022 -0.041 -0.722
exp2_data_contact_for_plot <- exp2_data_contact %>%
mutate(condition = case_when(
condition == "vh" ~ "Harmony",
condition == "vd" ~ "Disharmony"
))
exp2_data_contact_for_plot$condition <- factor(exp2_data_contact_for_plot$condition, levels = c("Harmony", "Disharmony"))
exp2_model_contact_for_plot <- glmer(data = exp2_data_contact_for_plot,
dominant_index ~ order * condition
+ (1 | participant) + (0 + order | participant) + (1 | testing_stem),
family = "binomial")
exp2_plot_contact <- plot_model(exp2_model_contact_for_plot,
type = "pred",
terms = c("order [all]", "condition"),
title = "",
line.size = 2.5) +
geom_hline(yintercept = 0.5,
linetype = "dashed",
size = 0.5,
color = "grey70") +
# the line colors
scale_color_manual(values = c(color_list[1], color_list[2])) +
# the error shade colors
scale_fill_manual(values = c(color_list[1], color_list[2])) +
scale_x_continuous(
limits = c(1, 48),
breaks = seq(0, 48, by = 8)
) +
labs(x = "Trial during contact",
y = "% Contact pattern",
color = "Contact\nlanguage") +
theme_bw(base_size = 25) +
theme(
axis.text.x = element_text(color = "black"),
axis.text.y = element_text(color = "black"),
legend.key.spacing.y = unit(0.5, "lines"),
legend.box.spacing = margin(0.5)
)
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
exp2_plot_contact
ggsave(plot = exp2_plot_contact,
filename = here("new_analysis_202406", "figures", "exp2_plot_contact.pdf"),
height = 6,
width = 8)
ggsave(plot = exp2_plot_contact,
filename = here("new_analysis_202406", "figures", "exp2_plot_contact.png"),
dpi = 300,
height = 6,
width = 8)
Setting the reference level: VH
exp2_data_testing$condition <- factor(exp2_data_testing$condition, levels = c("vd", "vh"))
exp2_data_testing$phase <- factor(exp2_data_testing$phase, levels = c("post_testing", "pre_testing"))
Setting up sum coding scheme
contrasts(exp2_data_testing$condition) <- contr.sum(2)
contrasts(exp2_data_testing$condition)
## [,1]
## vd 1
## vh -1
contrasts(exp2_data_testing$phase) <- contr.sum(2)
contrasts(exp2_data_testing$phase)
## [,1]
## post_testing 1
## pre_testing -1
exp2_model_testing <- glmer(data = exp2_data_testing,
dominant_index ~ condition * phase
+ (1 + phase | participant) + (1 | testing_stem),
family = "binomial")
summary(exp2_model_testing)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: dominant_index ~ condition * phase + (1 + phase | participant) +
## (1 | testing_stem)
## Data: exp2_data_testing
##
## AIC BIC logLik deviance df.resid
## 2538.4 2582.6 -1261.2 2522.4 1847
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5605 -1.0462 0.7145 0.8907 1.4888
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## participant (Intercept) 0.07346 0.2710
## phase1 0.03245 0.1801 0.59
## testing_stem (Intercept) 0.03533 0.1880
## Number of obs: 1855, groups: participant, 46; testing_stem, 42
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.14899 0.08071 1.846 0.0649 .
## condition1 -0.05537 0.07176 -0.772 0.4403
## phase1 0.10873 0.07497 1.450 0.1470
## condition1:phase1 -0.09505 0.06525 -1.457 0.1452
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cndtn1 phase1
## condition1 0.001
## phase1 -0.385 -0.009
## cndtn1:phs1 -0.009 -0.317 0.001
exp2_posthoc_testing <- emmeans(exp2_model_testing, pairwise ~ condition * phase, adjust = "mvt",
type = "response") %>%
summary(infer = c(T, T))
exp2_posthoc_testing
## $emmeans
## condition phase prob SE df asymp.LCL asymp.UCL null z.ratio
## vd post_testing 0.527 0.0292 Inf 0.469 0.583 0.5 0.915
## vh post_testing 0.601 0.0285 Inf 0.544 0.655 0.5 3.439
## vd pre_testing 0.520 0.0428 Inf 0.436 0.603 0.5 0.466
## vh pre_testing 0.500 0.0425 Inf 0.418 0.583 0.5 0.003
## p.value
## 0.3600
## 0.0006
## 0.6410
## 0.9973
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## Tests are performed on the logit scale
##
## $contrasts
## contrast odds.ratio SE df asymp.LCL asymp.UCL
## vd post_testing / vh post_testing 0.74 0.119 Inf 0.491 1.12
## vd post_testing / vd pre_testing 1.03 0.204 Inf 0.617 1.71
## vd post_testing / vh pre_testing 1.11 0.230 Inf 0.655 1.89
## vh post_testing / vd pre_testing 1.39 0.289 Inf 0.814 2.37
## vh post_testing / vh pre_testing 1.50 0.299 Inf 0.904 2.50
## vd pre_testing / vh pre_testing 1.08 0.241 Inf 0.612 1.91
## null z.ratio p.value
## 1 -1.875 0.2358
## 1 0.138 0.9991
## 1 0.516 0.9545
## 1 1.574 0.3895
## 1 2.051 0.1669
## 1 0.357 0.9842
##
## Confidence level used: 0.95
## Conf-level adjustment: mvt method for 6 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: mvt method for 6 tests
## Tests are performed on the log odds ratio scale
exp2_summary_testing <- emmeans(exp2_model_testing, pairwise ~ condition * phase,
type = "response")$emmeans %>%
tidy(conf.int = TRUE)
exp2_summary_testing
Setting the level order
exp2_summary_testing$condition <- factor(exp2_summary_testing$condition, levels = c("vh", "vd"))
exp2_summary_testing$phase <- factor(exp2_summary_testing$phase, levels = c("pre_testing", "post_testing"))
exp2_plot_testing <- ggplot(
data = exp2_summary_testing,
aes(
x = phase,
y = prob,
fill = condition
)
) +
geom_bar(stat = "identity",
color = "black",
position = "dodge",
width = 0.8,
alpha = 0.85) +
geom_errorbar(aes(ymin = conf.low,
ymax = conf.high),
position = position_dodge(0.8),
color = "black",
width = 0.2
) +
geom_hline(yintercept = 0.5,
linetype = "dashed",
size = 0.75,
color = "grey70") +
annotate("text",
x = 1.8,
y = 0.72,
label = "***",
size = 12) +
scale_fill_manual(values = c(color_list[1], color_list[2]),
labels = c("Harmony", "Disharmony")) +
scale_x_discrete(limits = c("pre_testing", "post_testing"),
labels = c("Pre-contact", "Post-contact")) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.25),
labels = scales::percent
) +
labs(x = "Testing phase",
y = "% Contact pattern",
fill = "Contact\nlanguage") +
theme_bw(base_size = 25) +
theme(
axis.text.x = element_text(colour = "black"),
axis.text.y = element_text(colour = "black"),
legend.key.spacing.y = unit(0.5, "lines"),
legend.box.spacing = margin(0.5)
)
exp2_plot_testing
### 3.9.3 Saving the plot
ggsave(plot = exp2_plot_testing,
filename = here("new_analysis_202406", "figures", "exp2_plot_testing.pdf"),
height = 6,
width = 8.5)
ggsave(plot = exp2_plot_testing,
filename = here("new_analysis_202406", "figures", "exp2_plot_testing.png"),
dpi = 300,
height = 6,
width = 8.5)
exp3_data_original <- read_csv(here("new_analysis_202406", "data_for_new_analysis", "cleaned_data", "exp3_new_cleaned_data_20240608.csv"))
## Rows: 4416 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (19): condition, chain, generation, participant, participant_id, phase, ...
## dbl (7): training_round, training_trial_in_round, training_trial_general, t...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
exp3_data_original
Two conditions; Four chains per conditions; four generations;
A total of 32 participants.
exp3_data_original %>% select(condition, chain, generation, participant_id) %>%
distinct()
exp3_unseen_only <- exp3_data_original %>%
filter(response_type == "unseen") %>%
mutate(
unseen_dominant_index = case_when(
condition == "vh" & unseen_type == "harmonic" ~ 1,
condition == "vd" & unseen_type == "disharmonic" ~ 1,
TRUE ~ 0
))
exp3_unseen_only
No difference between two conditions
exp3_unseen_only %>%
group_by(condition) %>%
summarize(total = length(unseen_type),
n_dominant_unseen = sum(unseen_dominant_index),
dominant_unseen_rate = n_dominant_unseen / total)
t.test(unseen_dominant_index ~ condition, data = exp3_unseen_only) %>%
tidy()
exp3_data <- exp3_data_original %>%
filter(phase == "testing") %>%
mutate(
dominant_index = case_when(
condition == "vh" & response_type == "harmonic" ~ 1,
condition == "vd" & response_type == "disharmonic" ~ 1,
TRUE ~ 0
)
)
exp3_data
Setting up the reference levels: VH for condition; gen1 for
generation
exp3_data$condition <- factor(exp3_data$condition, levels = c("vd", "vh"))
exp3_data$generation <- factor(exp3_data$generation, levels = c("gen2", "gen3", "gen4", "gen1"))
Setting up sum coding scheme
contrasts(exp3_data$condition) <- contr.sum(2)
contrasts(exp3_data$generation) <- contr.sum(4)
contrasts(exp3_data$condition)
## [,1]
## vd 1
## vh -1
contrasts(exp3_data$generation)
## [,1] [,2] [,3]
## gen2 1 0 0
## gen3 0 1 0
## gen4 0 0 1
## gen1 -1 -1 -1
exp3_model <- glmer(data = exp3_data,
dominant_index ~ condition * generation + (1 | participant_id),
family = binomial)
summary(exp3_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: dominant_index ~ condition * generation + (1 | participant_id)
## Data: exp3_data
##
## AIC BIC logLik deviance df.resid
## 1849.2 1896.0 -915.6 1831.2 1335
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3756 -0.9032 -0.7537 1.0853 1.4271
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.04027 0.2007
## Number of obs: 1344, groups: participant_id, 32
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.14259 0.06584 -2.166 0.0303 *
## condition1 -0.14356 0.06585 -2.180 0.0292 *
## generation1 -0.28544 0.11469 -2.489 0.0128 *
## generation2 -0.05012 0.11362 -0.441 0.6592
## generation3 0.06958 0.11371 0.612 0.5406
## condition1:generation1 0.05484 0.11466 0.478 0.6325
## condition1:generation2 0.14339 0.11363 1.262 0.2070
## condition1:generation3 -0.04982 0.11371 -0.438 0.6613
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cndtn1 gnrtn1 gnrtn2 gnrtn3 cnd1:1 cnd1:2
## condition1 0.002
## generation1 0.010 0.008
## generation2 -0.006 -0.001 -0.335
## generation3 -0.005 0.002 -0.335 -0.329
## cndtn1:gnr1 0.007 0.010 0.010 -0.004 -0.006
## cndtn1:gnr2 -0.001 -0.006 -0.005 0.001 -0.001 -0.335
## cndtn1:gnr3 0.002 -0.005 -0.006 -0.001 0.004 -0.335 -0.329
exp3_posthoc <- emmeans(exp3_model, pairwise ~ generation | condition, adjust = "mvt",
type = "response") %>%
summary(infer = c(T, T))
exp3_posthoc
## $emmeans
## condition = vd:
## generation prob SE df asymp.LCL asymp.UCL null z.ratio p.value
## gen2 0.374 0.0443 Inf 0.292 0.464 0.5 -2.732 0.0063
## gen3 0.452 0.0459 Inf 0.364 0.542 0.5 -1.041 0.2978
## gen4 0.434 0.0457 Inf 0.347 0.524 0.5 -1.433 0.1519
## gen1 0.458 0.0460 Inf 0.370 0.548 0.5 -0.911 0.3625
##
## condition = vh:
## generation prob SE df asymp.LCL asymp.UCL null z.ratio p.value
## gen2 0.416 0.0453 Inf 0.331 0.507 0.5 -1.820 0.0688
## gen3 0.452 0.0459 Inf 0.365 0.542 0.5 -1.040 0.2984
## gen4 0.530 0.0461 Inf 0.440 0.618 0.5 0.651 0.5152
## gen1 0.602 0.0449 Inf 0.512 0.686 0.5 2.214 0.0268
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## Tests are performed on the logit scale
##
## $contrasts
## condition = vd:
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## gen2 / gen3 0.723 0.192 Inf 0.366 1.428 1 -1.223 0.6119
## gen2 / gen4 0.779 0.206 Inf 0.394 1.538 1 -0.944 0.7810
## gen2 / gen1 0.706 0.187 Inf 0.358 1.393 1 -1.316 0.5528
## gen3 / gen4 1.076 0.283 Inf 0.548 2.112 1 0.280 0.9923
## gen3 / gen1 0.976 0.256 Inf 0.498 1.912 1 -0.093 0.9997
## gen4 / gen1 0.907 0.238 Inf 0.462 1.779 1 -0.373 0.9823
##
## condition = vh:
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## gen2 / gen3 0.863 0.227 Inf 0.439 1.697 1 -0.559 0.9442
## gen2 / gen4 0.631 0.166 Inf 0.321 1.241 1 -1.750 0.2977
## gen2 / gen1 0.470 0.124 Inf 0.238 0.928 1 -2.853 0.0222
## gen3 / gen4 0.731 0.191 Inf 0.373 1.433 1 -1.196 0.6296
## gen3 / gen1 0.545 0.144 Inf 0.276 1.072 1 -2.306 0.0963
## gen4 / gen1 0.745 0.196 Inf 0.378 1.466 1 -1.120 0.6773
##
## Confidence level used: 0.95
## Conf-level adjustment: mvt method for 6 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: mvt method for 6 tests
## Tests are performed on the log odds ratio scale
The VH dominant condition
exp3_data_gen1_vh <- exp3_data %>% filter(generation == "gen1" & condition == "vh")
t.test(exp3_data_gen1_vh$dominant_index, mu = 0.72) %>% tidy()
The VD dominant condition
exp3_data_gen1_vd <- exp3_data %>% filter(generation == "gen1" & condition == "vd")
t.test(exp3_data_gen1_vd$dominant_index, mu = 0.72) %>% tidy()
From gen1 to gen4
exp3_summary_by_chain_gen1_to_gen4 <-
exp3_data %>%
group_by(generation, chain, condition, participant_id) %>%
summarize(
dominant_count = sum(dominant_index),
dominant_rate = dominant_count / 42
)
## `summarise()` has grouped output by 'generation', 'chain', 'condition'. You can
## override using the `.groups` argument.
exp3_summary_by_chain_gen1_to_gen4
Adding the initial input
exp3_summary_by_chain_initial <- exp3_summary_by_chain_gen1_to_gen4 %>%
filter(generation == "gen1") %>%
mutate(generation = "gen0",
dominant_count = 23,
dominant_rate = 0.72,
participant_id = as.character(str_glue("{generation}_{chain}")))
exp3_summary_by_chain_initial
Combining the initial and subsequent generations
exp3_summary_by_chain <- rbind(exp3_summary_by_chain_initial, exp3_summary_by_chain_gen1_to_gen4)
exp3_summary_by_chain
From gen1 to gen4
exp3_summary_est_gen1_to_gen4 <- emmeans(exp3_model, pairwise ~ generation | condition, adjust = "mvt", type = "response")$emmeans %>%
tidy(conf.int = TRUE)
exp3_summary_est_gen1_to_gen4
Adding the initial input
exp3_summary_est_initial <- exp3_summary_est_gen1_to_gen4 %>%
filter(generation == "gen1") %>%
mutate(generation = "gen0",
prob = 0.72,
std.error = NA,
conf.low = 0,
conf.high = 0,
statistic = NA,
p.value = NA)
exp3_summary_est_initial
Combining the initial and subsequent generations
exp3_summary_est <- rbind(exp3_summary_est_gen1_to_gen4, exp3_summary_est_initial) %>%
rename(dominant_rate = prob)
exp3_summary_est
Setting the level order
exp3_summary_by_chain$condition <- factor(exp3_summary_by_chain$condition, levels = c("vh", "vd"))
exp3_summary_by_chain$generation <- factor(exp3_summary_by_chain$generation, levels = c("gen0", "gen1", "gen2", "gen3", "gen4"))
exp3_summary_est$condition <- factor(exp3_summary_est$condition, levels = c("vh", "vd"))
exp3_summary_est$generation <- factor(exp3_summary_est$generation, levels = c("gen0", "gen1", "gen2", "gen3", "gen4"))
Preparing the significance bars
exp3_signif <- data.frame(
condition = c("vh", "vh", "vd"),
start = c("gen0", "gen1", "gen0"),
end = c("gen1", "gen2", "gen1"),
y = c(0.82, 0.75, 0.8),
label = c("**", "*", "***")
)
exp3_signif$condition <- factor(exp3_signif$condition, levels = c("vh", "vd"))
exp3_plot <- ggplot() +
# individual chains
geom_point(
data = exp3_summary_by_chain,
aes(
x = generation,
y = dominant_rate,
group = chain
),
color = "grey80"
) +
geom_line(
data = exp3_summary_by_chain,
aes(
x = generation,
y = dominant_rate,
group = chain
),
color = "grey80"
) +
# the estimates
geom_errorbar(
data = exp3_summary_est,
aes(
x = generation,
ymax = conf.high,
ymin = conf.low,
group = condition
),
width = 0.5
) +
geom_point(
data = exp3_summary_est,
aes(
x = generation,
y = dominant_rate,
group = condition,
color = condition
),
size = 4
) +
geom_line(
data = exp3_summary_est,
aes(
x = generation,
y = dominant_rate,
group = condition,
color = condition
),
linewidth = 2.5
) +
scale_color_manual(values = c(color_list[1], color_list[2])) +
scale_x_discrete(limits = c("gen0", "gen1", "gen2", "gen3", "gen4"),
labels = c("Initial", "1", "2", "3", "4")) +
scale_y_continuous(
limits = c(0.2, 0.9),
breaks = seq(0.2, 0.9, by = 0.1),
labels = scales::percent) +
labs(x = "Generation",
y = "% Dominant pattern") +
geom_hline(yintercept = 0.5,
linetype = "dashed",
size = 0.75,
color = "grey70") +
geom_hline(yintercept = 0.72,
linetype = "dashed",
size = 0.75,
color = "grey70") +
geom_signif(
data = exp3_signif,
aes(
xmin = start,
xmax = end,
annotations = label,
y_position = y
),
textsize = 12,
manual = TRUE
) +
facet_grid(. ~ condition,
labeller = labeller(
condition = c(vh = "Harmony-dominant", vd = "Disharmony-dominant")
)
) +
theme_bw(base_size = 25) +
theme(
legend.position = "none",
axis.text.x = element_text(colour = "black"),
axis.text.y = element_text(colour = "black")
)
## Warning in geom_signif(data = exp3_signif, aes(xmin = start, xmax = end, :
## Ignoring unknown aesthetics: xmin, xmax, annotations, and y_position
exp3_plot
ggsave(plot = exp3_plot,
filename = here("new_analysis_202406", "figures", "exp3_plot.pdf"),
height = 6,
width = 10)
ggsave(plot = exp3_plot,
filename = here("new_analysis_202406", "figures", "exp3_plot.png"),
dpi = 300,
height = 6,
width = 10)